import string
import scipy
import Tkinter, tkFileDialog
import numpy as np
import pandas as pd
from scipy.spatial.distance import cdist
from scipy import ndimage
import skimage
import skimage.filters
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import matplotlib.image as mpimg
from PIL import Image
import cmath
import os
import sys
import glob
sys.path.append(os.path.abspath("C:\Users\Scherer Lab E\Documents\GitHub\Python_Data_Analysis"))
from common_functions import *
%matplotlib inline

Load the Wavefront Correction

wfc = mpimg.imread("K:\Pat's_Projects\Hyperuniformity\Exp08121501\90deg Shift non-ideal\\non-ideal_LUT_phase_correction.bmp")

Loading the Hologram and aperture

cd "J:\Pat's Projects\Dynamical Phase Transition\Phase Mask Scripts\Pattern\LSM\control_software\LCOS_Ctrl\data"

J:\Pat's Projects\Dynamical Phase Transition\Phase Mask Scripts\Pattern\LSM\control_software\LCOS_Ctrl\data

holo = mpimg.imread('Rn-L0R1.5.bmp')
holo = holo.astype(np.float32)
# Convert to Float
holo *= 2*np.pi/float(np.max(holo))
holo = holo[:, 96:-96]
plt.imshow(holo)
apeture = np.ones(holo.shape)
index_arr = np.indices(apeture.shape)
# Hologram Shape is 600 tall x 792 wide for Hamamatsu SLM
index_arr[0,:,:] -= 300
index_arr[1,:,:] -= 300
#index_arr[1,:,:] -= 792/2
# Draw the largest circle that will fit in the hologram
apeture[index_arr[0,:,:]**2 + index_arr[1,:,:]**2 >= (len(apeture[:,1])/10)**2] = 0

complex_arr = apeture*(np.cos(holo) + 1j*np.sin(holo))

''' If you want to look at either of the complex components'''
plt.imshow(np.real(complex_arr))
#plt.imshow(np.imag(complex_arr))

<matplotlib.image.AxesImage at 0x25e60198>

Load the CCD Data

cd "K:\Pat's_Projects\Hyperuniformity\Exp08181501"

K:\Pat's_Projects\Hyperuniformity\Exp08181501

ring_trap = mpimg.imread("Result of Pic_08181506 - Rn-1.5 ring trap_cropped_blur_thresh.tif")
# Resize data to roughly fit the size of the 4 times fft of the phase mask
ring_trap = ring_trap[:, 96:-96]
ring_trap = scipy.misc.imresize(ring_trap, 1/(5.0*1.8))
plt.imshow(ring_trap)

<matplotlib.image.AxesImage at 0x30139fd0>

Algorithm

def gs_phase_reconstruction(apeture, phase_mask, current_amp):
    '''
    :param apeture: 2D array of amplitude applied to phase mask
    :param phase_mask: 2D array of phases from 0 to 2*pi
    :param current_amp: The amplitude seen in the fourier plane that needs to
    be corrected for.
    '''
    init_complex_arr = apeture*(np.cos(phase_mask) + 1j*np.sin(phase_mask))
    height = len(init_complex_arr[:,0])
    width = len(init_complex_arr[0,:])
    diff_h = abs(height-len(apeture[:,0]))
    diff_w = abs(width-len(apeture[0,:]))
    pad_apeture = np.pad(apeture, ((diff_h/2,diff_h/2),(diff_w/2,diff_w/2)),
                             'constant', constant_values=0)
    
    # First FFT
    image_plane = np.fft.fft2(complex_arr)
    image_plane_amp = np.fft.fftshift(np.abs(image_plane))
    image_plane_phase = np.angle(image_plane)
    
    diff_h = abs(height-len(current_amp[:,0]))
    diff_w = abs(width-len(current_amp[0,:]))
    #print diff_h, diff_w, current_amp.shape
    pad_current_amp = np.pad(current_amp, ((diff_h/2,diff_h/2),(diff_w/2,diff_w/2)),
                             'constant', constant_values=0)
    pad_current_amp = np.sqrt(pad_current_amp)
    
    plt.figure()
    plt.subplot(221)
    plt.title('Ideal Amplitude')
    plt.imshow(image_plane_amp, cmap='gray')
    #plt.axis('equal')
    #plt.ylim(250,350)
    #plt.xlim(250,350)
    #plt.gca().invert_yaxis()
    plt.subplot(222)
    plt.title('Image Plane Phase')
    plt.imshow(image_plane_phase, cmap='gray')
    plt.subplot(223)
    plt.title('Current Amplitude')
    plt.imshow(pad_current_amp, cmap='gray')
    #plt.ylim(250,350)
    #plt.xlim(250,350)
    plt.gca().invert_yaxis()
    plt.subplot(224)
    plt.title('Ideal Phase')
    plt.imshow(holo, cmap='gray')
    plt.tight_layout()
    plt.show()
    
    for x in range(5):
        image_plane_amp = np.fft.fftshift(pad_current_amp)
        image_plane = image_plane_amp*(np.cos(image_plane_phase) + 1j*np.sin(image_plane_phase))
        
        slm_plane = np.fft.ifft2(image_plane)
        slm_plane_amp = np.fft.fftshift(np.abs(slm_plane))
        slm_plane_phase = np.angle(slm_plane)
        slm_plane_amp = apeture
        
        slm_plane = slm_plane_amp*(np.cos(slm_plane_phase) + 1j*np.sin(slm_plane_phase))
        
        image_plane = np.fft.fft2(slm_plane)
        image_plane_amp = np.fft.fftshift(np.abs(image_plane))
        image_plane_phase = np.angle(image_plane)
        
        plt.figure()
        plt.subplot(221)
        plt.title('Image Plane Amp')
        plt.imshow(np.abs(image_plane_amp), cmap='gray')
        plt.ylim(250,350)
        plt.xlim(250,350)
        plt.gca().invert_yaxis()
        plt.subplot(222)
        plt.title('SLM Plane Amp')
        plt.imshow(slm_plane_amp, cmap='gray')
        plt.subplot(223)
        plt.title('Image Plane Phase')
        plt.imshow(image_plane_phase, cmap='gray')
        plt.subplot(224)
        plt.title('SLM Plane Phase')
        plt.imshow(slm_plane_phase, cmap='gray')
        plt.tight_layout()
        plt.show()
        
    return slm_plane_phase

phase_correction = gs_phase_reconstruction(apeture, holo, ring_trap)

Apply GS Phase to WFC

def wrap2value(array, lower_value, upper_value):
    """
    Wraps the phase around 2pi such that multiples of 2pi greater than 2pi will equal 2pi and multiples of
     2p less than 0 will equal 0. This is called automatically on the phase attribute
    """
    array[array % (upper_value) != lower_value] %= (upper_value)
    array[np.logical_and(array % (upper_value) == lower_value, array / (upper_value) > lower_value)] = upper_value
    array[np.logical_and(array % (upper_value) == lower_value, array / (upper_value) < lower_value)] = lower_value
    return array

phase_correction_scaled = phase_correction+abs(np.min(phase_correction))
phase_correction_scaled *= 255/np.max(phase_correction_scaled)

new_wfc = wfc - phase_correction_scaled
new_wfc = wrap2value(new_wfc, 0, 255)
plt.imshow(new_wfc, cmap='gray')
plt.colorbar()
new_wfc = Image.fromarray(new_wfc.astype('uint8'))
new_wfc.save('L1_vortex_wfc_with_gs.bmp', format='bmp')